{
 "metadata": {
  "name": "",
  "signature": "sha256:3ab2e25d1c46652f9327cb6954f8f799d10cb26f9b4f9b5feb46b7735c56ea33"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Nonnegative matrix factorization\n",
      "\n",
      "A derivative work by Judson Wilson, 6/2/2014.    \n",
      "Adapted from the CVX example of the same name, by Argyris Zymnis, Joelle Skaf and Stephen Boyd\n",
      "\n",
      "## Introduction\n",
      "\n",
      "We are given a matrix $A \\in \\mathbf{\\mbox{R}}^{m \\times n}$ and are interested in solving the problem:\n",
      "    \\begin{array}{ll}\n",
      "    \\mbox{minimize}   & \\| A - YX \\|_F \\\\\n",
      "    \\mbox{subject to} & Y \\succeq 0 \\\\\n",
      "                      & X \\succeq 0,\n",
      "    \\end{array}\n",
      "where $Y \\in \\mathbf{\\mbox{R}}^{m \\times k}$ and $X \\in \\mathbf{\\mbox{R}}^{k \\times n}$.\n",
      "\n",
      "This example generates a random matrix $A$ and obtains an\n",
      "*approximate* solution to the above problem by first generating\n",
      "a random initial guess for $Y$ and then alternatively minimizing\n",
      "over $X$ and $Y$ for a fixed number of iterations.\n",
      "\n",
      "## Generate problem data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import cvxpy as cvx\n",
      "import numpy as np\n",
      "\n",
      "# Ensure repeatably random problem data.\n",
      "np.random.seed(0)\n",
      "\n",
      "# Generate random data matrix A.\n",
      "m = 10\n",
      "n = 10\n",
      "k = 5\n",
      "A = np.random.rand(m, k).dot(np.random.rand(k, n))\n",
      "\n",
      "# Initialize Y randomly.\n",
      "Y_init = np.random.rand(m, k)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Perform alternating minimization"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Ensure same initial random Y, rather than generate new one\n",
      "# when executing this cell.\n",
      "Y = Y_init \n",
      "\n",
      "# Perform alternating minimization.\n",
      "MAX_ITERS = 30\n",
      "residual = np.zeros(MAX_ITERS)\n",
      "for iter_num in range(1, 1+MAX_ITERS):\n",
      "    # At the beginning of an iteration, X and Y are NumPy\n",
      "    # array types, NOT CVXPY variables.\n",
      "\n",
      "    # For odd iterations, treat Y constant, optimize over X.\n",
      "    if iter_num % 2 == 1:\n",
      "        X = cvx.Variable(k, n)\n",
      "        constraint = [X >= 0]\n",
      "    # For even iterations, treat X constant, optimize over Y.\n",
      "    else:\n",
      "        Y = cvx.Variable(m, k)\n",
      "        constraint = [Y >= 0]\n",
      "    \n",
      "    # Solve the problem.\n",
      "    obj = cvx.Minimize(cvx.norm(A - Y*X, 'fro'))\n",
      "    prob = cvx.Problem(obj, constraint)\n",
      "    prob.solve(solver=cvx.SCS)\n",
      "\n",
      "    if prob.status != cvx.OPTIMAL:\n",
      "        raise Exception(\"Solver did not converge!\")\n",
      "    \n",
      "    print 'Iteration {}, residual norm {}'.format(iter_num, prob.value)\n",
      "    residual[iter_num-1] = prob.value\n",
      "\n",
      "    # Convert variable to NumPy array constant for next iteration.\n",
      "    if iter_num % 2 == 1:\n",
      "        X = X.value\n",
      "    else:\n",
      "        Y = Y.value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Iteration 1, residual norm 2.76585686659\n",
        "Iteration 2, residual norm 0.577758799504\n",
        "Iteration 3, residual norm 0.46343315761\n",
        "Iteration 4, residual norm 0.300312085357\n",
        "Iteration 5, residual norm 0.172468695929\n",
        "Iteration 6, residual norm 0.117552622713"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 7, residual norm 0.0855259222075\n",
        "Iteration 8, residual norm 0.0660380454036\n",
        "Iteration 9, residual norm 0.0530018181734\n",
        "Iteration 10, residual norm 0.0442728793651\n",
        "Iteration 11, residual norm 0.0364005958705"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 12, residual norm 0.0308842140499\n",
        "Iteration 13, residual norm 0.0256059616668\n",
        "Iteration 14, residual norm 0.0226869576657\n",
        "Iteration 15, residual norm 0.0191546943234\n",
        "Iteration 16, residual norm 0.0166449632154"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 17, residual norm 0.0135201384604\n",
        "Iteration 18, residual norm 0.0119471133563\n",
        "Iteration 19, residual norm 0.0149438374084\n",
        "Iteration 20, residual norm 0.0138663023673"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 21, residual norm 0.00922230392493\n",
        "Iteration 22, residual norm 0.00857605731059\n",
        "Iteration 23, residual norm 0.0074862441594"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 24, residual norm 0.00739813239648\n",
        "Iteration 25, residual norm 0.0100134191882\n",
        "Iteration 26, residual norm 0.00944406772568"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 27, residual norm 0.008678201611\n",
        "Iteration 28, residual norm 0.00873112072225\n",
        "Iteration 29, residual norm 0.00798267920957\n",
        "Iteration 30, residual norm 0.00846182763828"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Output results"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#\n",
      "# Plot residuals.\n",
      "#\n",
      "\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "# Show plot inline in ipython.\n",
      "%matplotlib inline\n",
      "\n",
      "# Set plot properties.\n",
      "plt.rc('text', usetex=True)\n",
      "plt.rc('font', family='serif')\n",
      "font = {'family' : 'normal',\n",
      "        'weight' : 'normal',\n",
      "        'size'   : 16}\n",
      "plt.rc('font', **font)\n",
      "\n",
      "# Create the plot.\n",
      "plt.plot(residual)\n",
      "plt.xlabel('Iteration Number')\n",
      "plt.ylabel('Residual Norm')\n",
      "plt.show()\n",
      "\n",
      "#\n",
      "# Print results.\n",
      "#\n",
      "print 'Original matrix:'\n",
      "print A\n",
      "print 'Left factor Y:'\n",
      "print Y\n",
      "print 'Right factor X:'\n",
      "print X\n",
      "print 'Residual A - Y * X:'\n",
      "print A - Y * X\n",
      "print 'Residual after {} iterations: {}'.format(iter_num, prob.value)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEbCAYAAAD0yNLXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3WuMXPdd8PGvr8FuYq/XtCmtQPa6LUUPSPU6aVGAatVd\npwVahLp2ApXKgyAb55Ggr+rEeVRaC4kn6wQJwRt8AakEKlo7CaIvELHdPMvtoY0dJ4CqFuIbCr2l\n9TWlcS72Pi9+52TOzJ7ZmTNz5rrfj3Q0c86cOed/dnbOb/53kCRJkiRJkiRJkiRJkiRJkiRJ0iBY\n1oNzjgAzwGVgS7JtT4P37AbOAqPJ+qHOJE2S1M/21ayfJALKYvt/ILM+C0yXnShJUv87DdyTWT+c\nLPVcrFmfBI6WnShJUv/bVLP+DPDJOvuOszCAjAM3Sk6TJGnAjAMnFnl9isixZI0RAWRdpxIlSWps\nZY/Oux64C9gJ3LvIfiNUKs5TaY5kFLhaftIkSYNisUr0KRYWYdXNgWzZsmUecHFxcXEpttSW9DRl\neStvatNIzfqBZMlzMWf/dH1B7uPMmTPMz88P7fKZz3ym52nw+rw2r2/4FipdKgrpdgBJcxTZ3EPa\nFyWvTuMU0V8kaxQ4Vn7SJElFdDuAnCByG9ncw3bgSGbbGNX9PA7WrE9RP8ciSeqSbleiXyECwu5k\nfSNR9vZgZp9JYAfweLK+J9l/mggup4EnupHYfjMxMdHrJHTUMF/fMF8beH1LVS+GMumk+aQ8T5LU\npGXLlkEL8aAXleiSpCFgAJEktcQAIklqiQFEktQSA4gkqSUGEElSSwwgkqSWGEAkSS0xgEiSWmIA\nkSS1xAAiSWqJAUSS1BIDiCSpJQYQSVJLDCCSpJYYQCRJLTGASJJaYgCRJLVk6ALIjRu9ToEkLQ1D\nF0CuXu11CiRpaRi6AHLpUq9TIElLw9AFkIsXe50CSVoaDCCSpJYMXQCxCEuSumPoAog5EEnqjpU9\nOu/u5PF24ATwyCL77gDGgCPAJWAGeAw4l7ezORBJ6o5eBJBZYE9m/WTyWC+IjCbvmQUuA/dQJ3iA\nORBJ6pZuF2GtBy7UbDsAPLjIe+aBESIXMgo8sdgJzIFIUnd0O4BsBPYBmzLbLhEBYjFXgfPNnMAc\niCR1R7eLsM4C41QHg+3AsQbvmwHS0DDGInUm5kAkqTt6UQfyXOb5CLCTCCr1HKe6zmM/EVAO5e1s\nDkSSuqNXrbBSh4EPsHjxVG2F+TGiGCw3gJw/v5e9e+P5xMQEExMTbSZRkobL3Nwcc3NzbR9nWftJ\nadkscBR4apF9RoiiqxGiHgRgKnlfXv3N/JveNM/3v19mMiVpuC1btgxaiAe96kg4TXXw2Fpnv3ng\nYSrBA6IO5Ey9A7/yCrz6ahlJlCQtphcBZIpojvsMlea5d2deHyMCDMAVFjb73QE8UO/gGzZYkS5J\n3dDtOpARIucB0f8jdSTzfJIIEo8n6weJnuuXgS1EJXrdviAbNkRF+q23lpVkSVKebgeQyzTO9Ryi\nuoL8CosPdVJldNQciCR1w9ANpjg6alNeSeqGoQsg1oFIUncMXQAxByJJ3TF0AcQciCR1x9AFEHMg\nktQdQxdA0ma8kqTOGroAYjNeSeqOoQsg5kAkqTuGLoCYA5Gk7hi6AGIORJK6o5fDuXfC/LVr89x8\nc4zIu2zYrk6SOmDQhnPvmJtuisU5QSSps4YugICdCSWpG4YygNiZUJI6r8wA8lCJx2qLFemS1Hmt\nzAcyScwqmLUMmAEebDtFJbApryR1XtEAshvYV+e1+TbTUhpzIJLUeUWLsB4k5iNfnrM8VW7SWmcO\nRJI6r2gAuUj1XOZZ97eZltKYA5GkzisaQB4m5itfl/PabPvJKYc5EEnqvKJ1IF8gAsUl4DKRIwHY\nCKwvMV1tMQciSZ1XNIAcAUaAx3Nem24/OeUwByJJnVc0gEwB24Bnc17b335yymFHQknqvKJ1IM8S\nxVd5jrSZltI4lIkkdV7RHMj/AY4RfUFOZrYvI+pGbi8pXW0xByJJnVd0+N4bi7w2D6xo8ji7k8fb\ngRPAI03sfxYYTdYP1UvD/Pw8N27A6tVw7RqsbKWvvSQtIa0O51709noFuCd5zBqh+Wa8s8CezHqa\nk6kXRPYBT1LpqDhLVNjnVeQDsHw5rF8Ply/DD/9wk6mSJBVSNIAcYJEbdxPWAxdyjrmP+gFkhuj9\nnjqWrC+ajrQprwFEkjqjaCX6BeB54D05rz3WxPs3EsFiU2bbJSIHk2c8Z9slFg7muIBNeSWps4oG\nkPuIcrJWJ4s9SwSF85lt24lcRZ5RKp0VU5eTx7ze8G+wM6EkdVbRAHKAuOHn9QOZafIYz2WejwA7\ngV119h2hUnGeSsNC7fYq5kAkqbNa6Ym+C9hM5BpOEjmCZcn2eq2j6jkMfIDqHEnW5ZxtaeDIzV/s\n3bsXgNOn4StfmeBjH5somCRJGm5zc3PMzc21fZxeNeOFaE11lMWHgR8ngtTyBtveSMP8fExL8qlP\nwU03we/8ToEUSdIS1K1mvBBFTnknKjIa7zTVwWMr+cVip1iYCxmlfp1JZadReOGFAimSJBVSNIA8\nQv3msxuaPMYUEQSOU6njuJtKABkjAkp6noNU9/uYov6cJJXEbIB/+ZcmUyRJKqxoJXraH2Md8NFk\nSVtDNVP/MULkPA4QzXEvAqeJOpXUJHBvZn0PEVSmiR7pp4EnGp3ISnRJ6qxWirD2U32DhwgI/6uJ\n916mcdA6xMJg1GiokwVsxitJnVU0BzIL3EbkRO5MlvuA9wIPlZu09pgDkaTOKlrrfpIIIEVf65Y3\nWmF94xtw223wrW/1OEWS1OdabYVVNAeyWKFQXxUYpTmQJJ5IkkpWtA7kWWJk3APEsCQAW4hOhKdK\nTFfb1qyJx5dfhrVre5sWSRpGrYxpdYSF858/BtzVfnLa9kYRFsDb3gYnTsDb397DFElSn+t2R8Ix\nokd42qnvXAvH6bh0ZkIDiCSVr1EAeQh4MGf7WSpFWH3LudElqXMaVaI/QINh0/uZc6NLUuc0U4R1\nmcht5I2MCzGIYur2tlNUIjsTSlLnNBNAdgFfYmEz3VFiOPZtRHDZWW7S2mdnQknqnEYB5Fnyx7ia\nJFpjjRCDIu4ErpSbtPaZA5GkzmlUB7ItZ9ss0fJqhMid3EkfBg8wByJJnVSkGe9mItcxTtSJ7CR/\nDo++YQ5Ekjqn2aFMpoEzRPA4CLyDPg8eYA5EkjqpUQBZT1SUH0nWdxKj7+Z5sqxElcUciCR1TqMi\nrHNEXccpInjk9ThfB2xPlr5iDkSSOqdRABkhmuiOUH8e8tHk9b4b99aOhJLUOY0Gz7pIBIhGtgEn\nKD48fNmqBlO8fh1Wr4bXXoPlvU6ZJPWpTs0H0sw85wDPFNi3a1asgFtugSt92chYkgZbK8O597Oq\nHAjA5s1w/Dhs2dKjFElSn+vWjIQDx4p0SeqMoQ8gNuWVpM4Y+gBiDkSSOqPMAPLREo9VGnMgktQZ\njfqBvKfJ42wE9gFPtJec8pkDkaTOaBRAThU4VpGOhDuA24A9Tew3RgylcgmYAR6jwBzso6Pwne8U\nSJkkqSmNAkg6UVQzzbv2N7HPJDEg43ZicMZGRonh42eTtNxDgeABUYT19a8XeYckqRmNAsgRYjbC\nZjzQxD5fSpaNxPAnjcwn+40C55tMRxWHM5GkzmgUQHYVOFa9OdPbdTVZWmIluiR1RpEJpbI2UZ2D\n2ELME7Kx3QTlmKEyH/sY8EiRN1uJLkmdUTSATFOZG6TW8TbTUu+Y2TqP/URAaXrcLXMgktQZRQPI\ng8DDwEmiBdVDRK7jfqLFVNlqK8yPEc2F6waQvXv3vvF8YmKC226bMAciSRlzc3PMzc21fZyig2ed\nJqazhcgNzBKV22NEMLm7yePMEkVg9WY3JHn9YvKY1oFMAUep3wFywWCK8/MxpPv3vw833dRk6iRp\nCenWYIpnM8+PEzf0dPu2oidvYJ7I7WQr0MdorvnvG5Ytsx5EkjqhaBHWFSIHMA/cRcwDAtEpsJmJ\np1L1It0YsBV4PDnXhZrXd9Bcc+EqaT3IW99a9J2SpHqKBpAZohL9MnGD3wMcTl57uIn3byVyLdPA\nBiI3cRx4Nnl9kggSjyfrB4Hdyfm2EMVmhYdLMQciSeUrY0KpMaKeosiwJ52yoA4E4MMfhl274CMf\n6UGKJKnP9XJCqbNE8OjL0XjBpryS1AllDee+mT6cEz1lEZYkla9oHciNnG3z9Pnc6uZAJKl8rQxl\n8gCVgLGRGF13G9EPpC+NjsLzz/c6FZI0XIoGkGfJH4tqnGjW25fMgUhS+YrWgdTrLHiKCCJ9yToQ\nSSpf0RzIujrbtxPNefuSORBJKl/RALLYnB+Fe4h3izkQSSpfK5Xou6hudXWRSl+QvmQORJLKV7T5\n7T76OKdBnZ7or74Ka9fCa6/F4IqSpIpu9URfLHh8oejJu2X1alizBl56qdcpkaTh0agIa5boKNjI\nBjozoVRp0mKsdfWaAUiSCmkUQO7vSiq6IK1I37Sp1ymRpOHQqAjrbLJPutwJfInoDzKabHsHMQ7W\nzs4ls31WpEtSuRrlQA7UrP8x8M6abWeJlllPUpnHo+/YlFeSytUoB1I7bMlGYFPOfuuJCZ/6ljkQ\nSSpX0X4gh4gcxwFiNsHLRBHWDmJmwb5lDkSSylU0gDxAzD64q2b7Y8B9paSoQ8yBSFK5WplQahdR\ngX4bUak+Sh+PxJsyByJJ5WplKBOIoqvaoUs+CjzRXnI6Z3TUHIgklalRAHlP8vhc8jhJfsfCZcQw\nJ30bQCzCkqRyNQog/xf4HpWmu8cW2beZHus9YxGWJJWrUQCZpHoI98tEh8HaW/EosL/EdJXOHIgk\nlatRAKmt5zhM9ETPs6f95HSOORBJKlfR4Xu3EkVVl4HzwAwxG+HTwO+XmrLW5A7nDnDjRozK+/LL\nsGpVl1MlSX2sW8O5HyJyJTuB3USHwovJcf64wHF2ECP9NmM3ME0Eq5kC56iyfDmMjMDlxeZUlCQ1\nrWgz3rNE8DgHnCRyImkHwpNNvH8SGCdyLWea2H8fMcbWU8n6LBFMWhpzK60HefObW3m3JCmraA5k\nhAgeEIHgcOa1Zqqov0SMr3WK5rJLM1SCB0QrsNpe8E2zHkSSylM0BzIKrCNyEFBp1jtN+XOij+ds\nuwRMtXpAW2JJUnmKBpA9RKAYIwZPPE4EkUkW7yPSilEW5mrSGox1wNXCB7Q3uiSVpmgAOU6MvpuV\nzpNedkfCESKIZKW3/1FaDCAWYUlSOVodC2sTkQt5iihWOrfo3q3Jay+VBpSW8hEWYUlSeVoJIEeJ\neogzxBAn24AjwAdoIVewiItELiQrXa97nr17977xfGJigomJiTfWR0fhP/+ztPRJ0kCam5tjbm6u\n7eMU7Tiyn8gFHCbqQ25Lto8DDwEfbPI4s0QwaDSHyEWqi7GmiH4h9c5TtyMhwGc/C089BY8+2mQq\nJWkJ6FZHwjFi7o/HqC5GOgXcXuA49RI6RrToSh2sWZ9i4TztTbMORJLK00oz3jxbm3z/ViIITAMb\niGKw48CzyeuTRC/1tKPgHio90ceA07QxZLx1IJJUnqIB5EvACSIXMErUe2wjeowfbOL9zybLI3Ve\nP5QsWfX2LcwciCSVp3CZFxE8asekOk5Mb9tri9aBfPObMD4O3/52F1MkSX2u1TqQVgIIRHHSOJEL\nOUGlCKrXFg0gL78cAypeuwbLWr1ySRoy3Q4gWZNEECmzCW+rFg0gAGvWwPe+B296U5dSJEl9rlut\nsPJsIVplPVTCsTrO4UwkqRzNBJDDwPVkuSfn9XRY9vvLSlQnWZEuSeVo1AprN9Gs9gqwnqhAP05k\ndR4A7s3sW/ZovB1hU15JKkejALKLmEAq7Zexj+g9fi/Rk/xZYpKppymxuW0nmQORpHI0CiCXqZ79\n7wHgBvAw0clv4JgDkaRytFKJfpn84PHJNtPSFeZAJKkcjQJI7Wi4EEVWef53m2npCnMgklSORkVY\nY0Trq+zcHBuACzX75QWavjQ6Ct/4Rq9TIUmDr5kirKvJfulypWZ9Of3RibAp5kAkqRzNVKLXG4G3\n1uk209IVdiSUpHI0yoHMFjhWM6Px9pyV6JJUjmEbUrDhWFinT8Odd8LZek0BJGmJ6eVYWAPFHIgk\nlWPJ5UCuX4fVq+HVV2HFii6lSpL6mDmQJq1YAevWwZUrvU6JJA22JRdAwKa8klSGJRlArAeRpPYt\nyQBiDkSS2rckA4g5EElq35IMIOZAJKl9SzKAOJyJJLVvyQYQi7AkqT1LMoBs2AAvvtjrVEjSYOtV\nANkNTAMzybKYHcD9wGZi3pHdyfOWvf/98NRT8OlPR890SVJxvQgg+4BniLnWDwFbiGBSzygxKvAZ\nYjbEM8C5dhLwznfCM8/AP/wDfPCD5kYkqRW9CCAzwFOZ9WPArkX2nydyHmNEMHmijES89a1w7Bi8\n730wPg7/+I9lHFWSlo5uB5DxnG2XgKkG77sKnC87MStXwu/9Hhw8CNPT8Mgj0GAsRklSotsBZBSo\nbUCbzre+bpH3zRDFXNNEHUipfuEX4MQJePxx+OVftoWWJDWj2wFkhIVT5KYBpd7UuceJupLHk2UL\njSveC/uxH4O//3vYtAm2bYs6EklSfY3mRC/b5ZxtaeCo17WvtsL8GFERfyhv5717977xfGJigomJ\niaYTt3o1/OEfws/+LHzoQ/C7vwv33QfLhm3WFElL2tzcHHNzc20fp9u3xnHgJNU5n7xtqREisIwQ\n9SAQ9SVH6+zfcEKpZv3Hf8COHfBTPwUHDsDNN5dyWEnqO4MyodQpFuZCRolcRZ554GEqwQOiNdaZ\n8pNW7V3vgi9/GW66Ce64w6FPJKlWL5rxHqS638cUcCCzPpZ5/Qpwoeb9O4AHOpa6jLVr4U//NPqK\n/OIvwn//dzfOKkmDoVel+7uJToFjRDPeP8m8NkMEiQ8m6+uBe4mcyxbgaer3BSmtCKv6oPAbvwHf\n+hZ88YtRVyJJw6LVIqxhqx7uSAABeP316Cuydi187nOwfEmOIiZpGA1KHcjAWrkSPv95+OY34ROf\nsMOhJBlAClizJoqw/umfoomvJC1l3e4HMvDWr4e//dvoK7JxI/zWb/U6RZLUGwaQFtx6Kxw9Cj/3\ncxFEfvVXe50iSeo+A0iLNm+OnMjkZExQ9aEP9TpFktRd1oG04Sd/Ev7qr+DXfg3++Z97nRpJ6i4D\nSJvuuAP+7M9iFN+vfrXXqZGk7jGAlODnfx7+4A+iGOv8+V6nRpK6wzqQknzsY3DhAmzfDk8+CWNj\nvU6RJHWWAaREv/3b0UP9Z34mJqe6445ep0iSOsehTDrgb/4Gfv3X4Y/+CH7lV3qdGklanGNhhb4I\nIAD/+q/wkY/APffApz7lpFSS+pcBJPRNAIEYvfeXfgl+4ifg0KGYW0SS+o2DKfahH/kR+Lu/i3lE\ntm+PSnZJGhYGkA5buxaOHIkK9Z/+afj3f+91iiSpHAaQLli+HGZnYc8eeP/7oYS57CWp5wwgXfSb\nvwl/+Zdw993w2c/2OjWS1B4r0Xvg61+POdbvvhs+/Wn4oR/qdYokLWVWog+Qd78bvvzlaOr79rfD\nrl0xSdUAxD5JeoM5kB574YWYY/3RR+GVV+DjH49ly5Zep0zSUmE/kDBwASQ1Pw/PPAN//udRT/Ku\nd0UgueuumG9EkjrFABIGNoBkvfZaTFb16KMx8+Gdd0YwmZqKZsGSVCYDSBiKAJJ16VL0I/mLv4CT\nJ+HHfxze+1543/tiefe7YcWKXqdS0iAzgIShCyBZ167Bc8/B00/DV74Sy4svwm23RTBJA8vb3tbr\nlEoaJAaQMNQBJM+FC5WAkj6uXRvT7Y6Nxdzt6TI2BiMjvU6xpH4zaAFkN3AWGE3WD5W0/5ILILXm\n5+HMGfja1+Dcucpy9mw8rlxZHVA2b4Yf/VG49VZ4y1tisZ5FWloGKYDsA54EnkrWZ4ETwOMl7L/k\nA8hi5ucjx5INKOfORVPiF1+sLKtWVYLJW95SHVze/GZYvz6Wdesqz2+5JYZskTR4BimAXKSSkwCY\nBB4A7ixh/6EOIHNzc0xMTHT0HPPz8NJL8J3vVAJK9vl3vwtXrsDVq/GYLj/4Adx8cyWgZAPLunXx\nmH1e+3jLLfDcc3NMTU2wZk0EsWHSjc+ul7y+wdZqAOn2lLbjOdsuAVMl7T/UuvFPvGxZ3NTXrYN3\nvrP5912/Xh1U0ucvvRTP08cXXqis17723e/OsXz5BD/4QeRm1qyJ4rQ1a6qfZ7c1s9x0UwSk7LJy\n5cJtq1bB6tXVy6pV5UwGNuw3IK9vaep2ABklchRZl5PHdcDVNvdXj6xYER0e2+n0uHdvLPPz0Rfm\n5ZcjZ/Pyy9XPs9tql5deipxSdtsrr8TxapfXX1+47dVXK4/p81WrIgjlBZc0GKVL7Xq67Wtfi2LD\nNEhlH/Oer1gRQXTFisWfr1iRfx3Za8kuEAFx+fLKUrueblu1KsZpa2Z56aWYQC1972JLo/O2G7Dn\n5+HGjerPOO9xfr5xWtPlyhX4r/+qpLn2emr3LyLv71NvPS1gmZ+vXmq3Qfx46nSxcrcDyAjVxVFQ\nCRCjLAwIRffXEFi2rHKTXr++t2lJb0RpQKldXn+9sqQ3prxtn/88TE4uvLGnz69di1xYuu3GjcjV\nXb9e/bx2/caNCCJ5OanabWvWxN/2xo3KTfb116vXs8trr0XwvXat8XLhQgzJkx6n9gaXXbLnqz1v\nelNvJqDkbb9+PdK9fHl1TjMvwKc35GaWK1fg8OGF11BvvUgQaeZ42b9Neuy8oJVd/7d/i4Yyw2SK\nhTmKMeAGkaNod//TwLyLi4uLS6HlNC3odg7kIpGryErX83ITRfd/R+tJkyQV0e2Gl6eo1GGkRoFj\nJe0vSeqSXoyitAFYD3wtWb8X+GJmfQz4QGa90f6SpB7odU/0MaJZ7p9kXpsBdgAfbHJ/SZK0ROwg\nRhTIsxuYJn5IzHQtReWqd307gPuBzURd3u7kuTSQul2J3ilFx9YaJDuInNcRIvc1AzwGnOtlolo0\nSXQO3Q6cyXk9b9iaaeoPc9NvGl3fKHFNs0Td3j0M3ue4O3m8nRhS6JGc1wf5u7jY9Q36d3GESPNl\nIJ3zdE/NPoP++RW2j6gzSaU3nWFxL9Fs+QbRKu2jvU1OKWaB/Tnba5tsTwJHO5+c0tW7vhmi+fmm\nrqamPLW5qpNUbrgw+N/FRtc36N/FfTXrJ6nO5Q/659eSYbnp1DPoN508eTfYcRZ+luPEl3XQLBZA\nBtV6qm+mENeT/cwG+bvYzPUN+nfxNJHrTR1OllThz2/Qx09dKmNlXQXO9zoRHdZo2JphMUP8qptm\n4Q2rn20kfqFuymy7RKVf1qB/FxtdX2qQv4tTVDdA2gI8nTxv6fMb9DqQpTJWVvaX0BgLy52HwVIY\ntuY41eXl+4nPdhDKmc8SN5nzmW3bqfTJGvTvYqPrSw3yd/F85nmau//9ZH3QP7+W7GDhRY8Qf5hN\nXU9NZ9S20klvOoMsr4in6LA1/axeEVataVocQqIPjBCf16Zkfdi+i7XXB8PxXVxPpPkosDWzvaXP\nb9CLsGp7qUPlV2ztH2NQ1bbwOEbMhzJsig5bM2jSL2M2GF4hguQgOkxUuJ5P1oftu1h7fTAc38Ur\nRI73TqI4Kw2ALX1+gx5AvOkMj2EftmYeeJjq/8sx8pv79ru0KfJzmW3D9F3Mu75h+C7Wfj77gQPJ\n85Y+v0EPIN50BlO9ERAOUt1scIrKP/ggybu+K8CFmm07GLxfsNNE8UfaVyctBhmW72K96xv072Ja\nRJwNgOn/6Tpa/Px6MRZW2YZ5rKxXiH/g/5fZ9hDx62gQr28r8PFk2UZ8KV8Gvp28fhz4MPAjyeM3\ngc91P5kta3R9XwU+AfwP4C7gr4Enup/Mlk0BbyM+pzXJ8/+ZrMPgfxcXu75B/y5eID6f7P/bHuJ/\nMv2OFf78ejUWVtmGeays9cQHmfYefZrBuuloOKSVyrWOAHdn1gf1u9jM9Q36d3ErlWa5G4kfOA/W\n7DOon58kSZIkSZIkSZIkSZIkSZIkSZKkHpsiRp3Nzub2ZOb1yV4kqs55x4kOVO/pwvn3EX+L9G+y\nvub1Z6j8zW4QPZ877VhyrkEdJVjSkDoKXGfh8NBHWXjz7Ia8844Bz9OdAJJKg8jhnNdm6P6Q9oeJ\nv4Ek9Y0j5M8v0Ks5P/plrpHTxPzUN1g4B/UY3c8NHMEAoiYM+mi8GnzpaLvdHpetV+etZ2fyeITq\niYsukj9Xg9RzBhD10j4q9RBHiCKlyZrXTybb99dsT+sFNgP3E8Nqp3OMj2SOl/66z/6yn61z3v1U\nipNq60f2JfudTN6fup+oM0nfk9YhnGThDHaLOUd1EKkne+1p7ilbV5LdL03XdGaftNjuXuJvc5H8\nudmXEcVnJ5N9ZnP2afXzkaTC8oqw6pXxHwBOZNYPEzeq1HTyvlniJnWaSrHLsZp9Z6m+uS523slk\n+wdq0p2t9D9J9Q0zTcvF5H3rqdx4m5EtotpP5bogguHJJtK+m4XXmKbreeJvvjWTzi8k70/flw2Y\n6ef0yWR9c7KeveZ2Ph9JKiwvgNzLwpthOhf6PZltO2rem65/IVnfSqXi+zDVN/za99Y7L0QrrGwA\nSdPy0cw+6U03fW96/E9m9klv8tlz1lNbx5G2WJskP4Dkpf1+FgaQKRb+HY8QDRmy761t4ZXuk5UG\nNmj/89GQWNnrBEg5xpPH+4iJlyBmRztD3LzOZ/ZNZ0x7NrPtLiIXsIO4id6WbK+dsrNIWs5mtqXn\nmqJ6PohTmeeXWjhXaidR5HSEStrbURuALtN4mtmzdda3EnNhQOufj4aEAUT9Zpy4GUH8qn5qkX1h\n4c0RInDMZpZtLF6vkJ73VM720Zxt9dTedFv1LLCLKCbKa9qbZ2OB4zcqWhuhfuOCS7T/+WhIWImu\nXqu9mR2kUra+rea1HSysmK5toTRG3HRnidnUzpN/c807b550utYtmW1jNa+1a5SFRWmHkuOPEzPH\nZeUFgLw1WmnaAAABM0lEQVT92knPWM22Lcnxz1MJCq18PhoiBhB1U1qElP11m95gthNFQmeIX+DH\nidY8W5PXx4k5nM/VHLP2l3J6jrToZ4z4pQzVQSDvvFD5dZ0e9yzwGJEjSKW5g9pioA0slLcta1uS\n5u05r+0k/wac5nTS94wnx1lG9Q087+8N+QE1u88IkdOYyazvBB5I1k/R+ucjSYVkhzK5zsKhTA4n\nr52g+pd42qw23T99LW2Cep34NTxDtdnkPaeptDY6ClygunVV7Xl3ZI77PNUV57NUmvFmK5yz7zlB\nVHzP5GzLcyDzN7lR8zdJTVLd2im1OznH0SRtacX+BaIF1HROGtKK8OvJubYSdRTXk/el1/V88vfY\nn7nmbAOBVKufjyRJkiRJkiRJkiRJkiRJkiRJkiRJkvrd/wfy3fMwDC58cwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1043ac910>"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Original matrix:\n",
        "[[ 1.323426    1.11061189  1.69137835  1.20020115  1.13216889  0.5980743\n",
        "   1.64965406  0.340611    1.69871738  0.78278448]\n",
        " [ 1.73721109  1.40464204  1.90898877  1.60774132  1.53717253  0.62647405\n",
        "   1.76242265  0.41151492  1.8048194   1.20313124]\n",
        " [ 1.4071438   1.10269406  1.75323063  1.18928983  1.23428169  0.60364688\n",
        "   1.63792853  0.40855006  1.57257432  1.17227344]\n",
        " [ 1.3905141   1.33367163  1.07723947  1.67735654  1.33039096  0.42003169\n",
        "   1.22641711  0.21470465  1.47350799  0.84931787]\n",
        " [ 1.42153652  1.13598552  2.00816457  1.11463462  1.17914429  0.69942578\n",
        "   1.90353699  0.45664487  1.81023916  1.09668578]\n",
        " [ 1.60813803  1.23214532  1.73741086  1.3148874   1.27589039  0.40755835\n",
        "   1.31904948  0.3469129   1.34256526  0.76924618]\n",
        " [ 0.90607895  0.6632877   1.25412229  0.81696721  0.87218892  0.50032884\n",
        "   1.245879    0.25079329  1.25017792  0.72155621]\n",
        " [ 1.5691922   1.47359672  1.76518996  1.66268312  1.43746574  0.72486628\n",
        "   1.97409333  0.39239642  2.09234807  1.16325748]\n",
        " [ 1.18723548  1.00282008  1.41532595  1.03836298  0.90382914  0.38460446\n",
        "   1.213473    0.23641422  1.32784402  0.27179726]\n",
        " [ 0.75789915  0.75119989  0.99502166  0.65444815  0.56073096  0.341146\n",
        "   1.02555143  0.24273668  1.01035919  0.49427978]]\n",
        "Left factor Y:\n",
        "[[  7.38991833e-01   3.15957978e-01   8.46211348e-01   7.90522539e-01\n",
        "    8.82326030e-01]\n",
        " [  6.37868033e-01   8.22907024e-01   5.32198000e-01   5.70689637e-01\n",
        "    6.21191813e+00]\n",
        " [  5.59748656e-01   6.34112010e-01   7.99615283e-01   1.72054035e-01\n",
        "    6.92576630e+00]\n",
        " [  2.61288516e-01   9.41947419e-01   4.03583183e-02   1.09118729e+00\n",
        "    9.07778543e-07]\n",
        " [  7.89189550e-01   3.41453292e-01   1.17654458e+00   3.93009044e-01\n",
        "    5.50024762e+00]\n",
        " [  7.39615442e-01   4.74493175e-01  -2.23332571e-04   6.74749299e-01\n",
        "    8.42579458e+00]\n",
        " [  4.73914127e-01   3.70454244e-01   8.08948369e-01   1.36848129e-01\n",
        "    3.44366220e-06]\n",
        " [  5.88504809e-01   7.27646377e-01   1.00390505e+00   1.03542480e+00\n",
        "    3.71366168e-01]\n",
        " [  8.14822860e-01   8.87015769e-04   2.91164377e-01   1.17787451e+00\n",
        "    1.24901335e-06]\n",
        " [  4.22680617e-01   7.77641517e-02   5.87259008e-01   6.51086033e-01\n",
        "    2.66173216e+00]]\n",
        "Right factor X:\n",
        "[[  1.13055890e+00   4.05899679e-01   1.59181960e+00   6.82867774e-01\n",
        "    9.75411818e-01   3.23464160e-01   8.83710480e-01   1.64529269e-01\n",
        "    9.23391090e-01   1.03847861e-01]\n",
        " [  9.03465524e-01   6.86715676e-01   6.62169881e-01   1.12490745e+00\n",
        "    1.03933855e+00   3.06370001e-01   7.29180054e-01   1.18625225e-01\n",
        "    8.75435486e-01   8.00971786e-01]\n",
        " [  6.63783207e-03   1.79085385e-01   3.11550072e-01   2.62447584e-02\n",
        "    1.60660298e-02   2.77495461e-01   6.46185026e-01   1.51538848e-01\n",
        "    5.43725876e-01   4.58269799e-01]\n",
        " [  2.23736169e-01   5.25847565e-01   2.32705796e-02   4.01864284e-01\n",
        "    8.80884850e-02   3.23502266e-02   2.59210460e-01   4.94700824e-02\n",
        "    3.53704441e-01   4.54771630e-02]\n",
        " [  2.28693369e-02   2.99295139e-02   2.74091017e-02   5.18043712e-04\n",
        "    1.64409925e-04   2.93841883e-04   1.72556068e-02   1.61648660e-02\n",
        "    6.42191908e-04   3.35760370e-02]]\n",
        "Residual A - Y * X:\n",
        "[[ -1.68636983e-04   3.60207904e-05  -3.98077462e-04  -2.04516172e-04\n",
        "   -4.16031924e-04   1.58422455e-03  -7.36906450e-04  -5.91624286e-05\n",
        "   -5.46879309e-04  -4.00508769e-04]\n",
        " [ -6.83414765e-04  -6.95866828e-04  -6.36012200e-04  -5.76770312e-05\n",
        "   -1.33274006e-04   6.27824831e-05  -3.32923224e-04  -3.46132086e-04\n",
        "    2.01209611e-04  -6.50509595e-04]\n",
        " [ -7.73470728e-04  -9.20826799e-04  -6.31154497e-04   2.22522673e-05\n",
        "    9.78991659e-05  -1.17532606e-03   8.50142570e-05  -4.05047780e-04\n",
        "    5.07774407e-04  -5.64939614e-04]\n",
        " [ -3.10948903e-04  -2.61145596e-04  -3.80145546e-04  -2.41030729e-04\n",
        "   -2.44643018e-04   4.30377374e-04  -2.61672784e-04  -1.20679505e-04\n",
        "   -2.79346735e-04  -4.08818401e-04]\n",
        " [ -7.07080893e-04  -8.12121976e-04  -6.37850976e-04  -4.46396107e-05\n",
        "    4.76451453e-05  -1.27571378e-03   9.33258194e-05  -3.50044450e-04\n",
        "    3.29244878e-04  -4.88627114e-04]\n",
        " [ -3.85692817e-04  -8.61465250e-04  -6.94274594e-04   5.50232944e-04\n",
        "    4.82360881e-04  -1.29338044e-03  -6.97353986e-04  -6.10252918e-04\n",
        "    2.71414348e-04  -1.10507632e-03]\n",
        " [ -3.89135610e-04  -3.02813339e-04  -7.79701822e-04   3.94758648e-04\n",
        "   -1.51231975e-04   4.63196558e-03  -1.25457666e-03  -4.81723061e-04\n",
        "    1.10896155e-05  -1.32234306e-03]\n",
        " [ -3.69191741e-04  -3.38504643e-04  -4.69486964e-04  -3.62561200e-04\n",
        "   -2.38561015e-04  -6.07402720e-04  -6.92616183e-05  -1.03315580e-04\n",
        "   -4.01477226e-04  -2.98251672e-04]\n",
        " [ -2.37003094e-04  -5.11457693e-05  -4.34540394e-04  -3.83365150e-05\n",
        "   -3.15650926e-04   1.86541629e-03  -7.04768572e-04  -1.45497987e-04\n",
        "   -2.65728261e-04  -5.29103341e-04]\n",
        " [ -6.65085219e-04  -9.73758356e-04  -3.70136974e-04  -1.03883251e-04\n",
        "    3.94287162e-04  -4.20725740e-03   1.14398039e-03  -2.59815037e-04\n",
        "    6.72891426e-04  -4.62543141e-06]]\n",
        "Residual after 30 iterations: 0.00846182763828\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}